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Abstract. In this article we propose a discrete lattice model to simulate the elastic, 
plastic and failure behaviour of isotropic materials. Focus is given on the mathematical 
derivation of the lattice elements, nodes and edges, in the presence of plastic deforma¬ 
tions and damage, i.e. stiffness degradation. By using discrete calculus and introducing 
non-local potential for plasticity, a force-based approach, we provide a matrix formulation 
necessary for software implementation. The output is a non-linear system with allowance 
for elasticity, plasticity and damage in lattices. This is the key tool for explicit analysis 
of micro-crack generation and population growth in plastically deforming metals, leading 
to macroscopic degradation of their mechanical properties and fitness for service. An il¬ 
lustrative example, analysing a local region of a node, is given to demonstrate the system 
performance. 
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1 Introduction 

Lattice models for analysis of deformation and fracture of solids have been developed over 
the last thirty years mainly for quasi-brittle materials, such as concretes and rocks @1, i, 
where local behaviour is elastic-brittle, i.e. the only mechanism for energy dissipation is 
the generation of new surfaces (brittle micro-cracking). Lattice models contain set of sites 
connected by bonds; in the language of algebraic topology this is a 1-complex embedded 
in or i.e. 2D or 3D graphs. The challenge with lattices is that they intend 
to represent a continuous solid by a discrete system. Specifically, the stored energy in 
any lattice region is required to be equivalent to the stored energy in the corresponding 
continuum region. This is used to derive a link between properties of lattice elements, 
e.g. bond stiffness coefficients, and macroscopic properties of the material, e.g. elastic 
modulus and Poisson’s ratio. For irregular graphs this is not possible exactly; approximate 
calibrations of lattice properties for known macroscopic properties are used, e.g. [TB]. For 
graphs with some regularity, the link can be established rigorously. For example, isotropic 
materials can be represented exactly by 2D graphs based on hexagonal structure, e.g. 
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m. and by 3D graphs based on truncated octahedral structure m- Other possible 
regular 3D graphs can represent cubic elasticity, i.e. three parameter elastic behaviour, 
but not isotropic elasticity, i.e. two parameter behaviour, typical for most engineering 
materials El- Nevertheless, lattices are being used because of the need to represent 
failure in materials, albeit not always with exactly calibrated local properties. Failure is 
a generation of a new internal surface, i.e. a non-topological concept. Hence, the classical 
continuum mechanics does not work, as it is a thermodynamic bulk theory. 

Our work addresses cleavage fracture, which is a phenomenon of fast and catastrophic 
crack propagation. It can be observed in all metals with body-centred cubic crystal lattices 
due to the smaller number of slip planes, i.e. reduced ability to dissipate energy by 
plastic deformation. We are specihcally interested in ferritic steels, used for example to 
manufacture nuclear reactor pressure vessels, where cleavage is a potential cause for serious 
concern. Cleavage is known to be triggered by rupturing of brittle second-phase particles, 
such as carbides, which typically decorate the grain boundaries of the polycrystal and can 
rupture due to plastic overload from the surrounding grains. The existing modelling of this 
phenomenon, the so called local approach to cleavage fracture, is based on the weakest-link 
statistic, i.e. global failure probability derives from the failure probabilities of individual 
particles, following prescribed size distribution but treated as independent events [3]. This 
approach works well for predicting cleavage fracture at very low temperatures, where 
plasticity is limited, the number of micro-cracks formed prior to cleavage is small and 
they are spatially scattered, i.e. only a set of largest particles failed. The probability 
of cleavage fracture, however, decreases rapidly with the increase of temperature, where 
plasticity is enhanced. This apparent increase of material toughness with temperature 
cannot be captured with the existing modelling strategy, despite of the many improvements 
in the particle rapture criterion over the years 0, m- However, accurate assessment 
with reduced conservatism is needed for more economic exploitation of reactors, including 
better planning of inspection intervals and life-extension decisions. From the current status 
of modelling it can be deduced that the plasticity enhancement with temperature leads 
to generation of increasingly large micro-crack populations, in which case the argument 
behind the weakest-link statistic is violated, e.g. there is a significant effect of micro¬ 
crack interactions prior to ultimate cleavage. This will explain the reduced probability 
of cleavage, hence increased toughness, of the material with temperature. Lattice models 
are particularly suitable to study the generation and interaction of micro-cracks as they 
grow, coalesce and progress to final failure. To this end we present a lattice development 
relevant to plastically deforming metals. 

A mathematically rigorous treatment of lattices can be achieved when they are anal¬ 
ysed as graphs [5]. However, the discrete exterior calculus (DEC) specialized to graphs in 
this reference, is developed for scalar problems. This means that the nodal unknown (a 
0-cochain) is scalar, i.e. temperature, pressure, concentration, etc. The gradient of this is 
also a scalar field over the edges (1-cochain). Things go relatively simply for such physical 
problems. In mechanics this approach is rather more difficult as we have a vector-valued 
nodal field, namely a displacement vector assigned to each node. Article [15] is the only 
one to our knowledge which attempts to apply DEC to mechanical problems in elastic 
settings. Our aim is to build upon this theoretical basis and present a graph-theoretical 
approach to lattices including elastic, plastic and damage behaviour. 
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Microstructure representation 


Micro-structures of metals and alloys develop by a crystallization process, which starts 
at spatially randomly distributed nuclei and finishes when neighbouring growing crystals 
touch each other to form grain boundaries [1]. Under uniform temperature distribution, 
the individual crystals grow with equal rates and therefore the final micro-structure is de¬ 
scribed by the Voronoi diagram constructed around the set of nuclei [2] . The 3D Voronoi 
diagram is constructed by the intersection of planes bisecting normally the segments con¬ 
necting each pair of neighbouring nuclei, as illustrated in Figure 1. Thus, each crystal is 
a polyhedron representing the neighbourhood of a nucleus, where all points are closer to 
the nucleus than to any other nucleus. Initially, we seek to represent the poly-crystal with 




Figure 1: Illustration of Voronoi diagrams constructed around a set of points: 
2D construction (left) providing comprehensive picture of the process; 3D 
construction (right) with the polyhedron formed around the central point 
shown in green. 


a regular tessellation of the 3D space, shown in Figure 2, where a grain is topologically 
equivalent to the average grain in arbitrary Voronoi tessellation [9]. This simplifies the 
geometry and allows for exact calibration of the lattice emerging from the tessellation [ni¬ 
ce// complexes and their duals 

Irrespective of whether we use an irregular tessellation of space (Voronoi around nuclei) or 
the regular tessellation of Figure 2, the material is subdivided into compactly packed and 
non-intrusive 3D cells. Using terminology from algebraic topology, see [7], the polycrystal 
is a p-complex, where p = 3 embedded in where d = 3. Thus, each crystal is a 3-cell, 
bounded by a number of 2-cells (planar polygonal faces), each of which is bounded by a 
number of 1-cells (edges), each of which is bounded by two 0-cells (nodes). It should be 
noted that p = d is not a requirement for p-complexes embedded in p-complexes with 
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Figure 2: Truncated octahedron (left), representing an average single grain. 
Assembly of cells, filling compactly a 3D region (right), representing topo¬ 
logically averaged polycrystal. 


1 < p < d can be formed and treated similarly. For example, the collection of all faces of 
the regularized polycrystal (or grain boundaries of the real microstructure) without the 
cell (grain) interiors is a 2-complex embedded in It contains 0-cells, the vertices of 
the truncated octahedrons (or quadruple points in the real microstructure), 1-cells, the 
edges of the truncated octahedrons (or triple lines), and 2-cells, the faces of the truncated 
octahedrons (or grain boundaries). The structure of thin-walled closed cell foam can be 
viewed as a 2-complex. Similarly, the collection of all edges of the regularized polycrys¬ 
tal (or triple lines in the real microstructure) without the cell (grain) interiors and faces 
(boundaries) is a 1-complex embedded in It contains 0-cells, the vertices of the trun¬ 
cated octahedrons (or quadruple points in the real microstructure), and 1-cells, the edges 
of the truncated octahedrons (or triple lines). The structure of thin-ligament open cell 
foam can be viewed as a 1-complex. Illustration of the three possibilities for complexes 
embedded in is given in Figure 3. The historical name for 1-complex is a graph, in this 
particular case the graph of edges in the right figure. 

Important for the analysis on complexes is that each p-complex can be used to build a 
dual p-complex, where the s-cells of the dual correspond to (p — s)-cells of the original or 
primal, for each 0 < s < p. Specifically for a 3-complex, the 0-cells of the dual correspond 
to 1-cells of the primal, the 1-cells of the dual correspond to 2-cells of the primal, the 2-cells 
of the dual correspond to 1-cells of the primal, and the 3-cells of the dual correspond to 0- 
cells of the primal. For a 3-complex embedded in ZD, such as the material microstructure, 
this means that the nodes of the dual will correspond to the grains, the edges of the dual 
will correspond to grain boundaries, the faces of the dual will correspond to triple lines, 
and the cells of the dual will correspond to quadruple points. 

Considering the representation of the material as a 3-complex, the topology of the 
site-bond model is derived as a 1-complex (graph) by placing sites (graph nodes) at the 
centres of all 3-cells and bonds (graph edges) between neighbouring sites as illustrated 
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Figure 3: Illustration of 3-complex (left), 2-complex (centre), and 1-complex 
(right), embedded in Illustration based on the topologically averaged 
poly crystal. 


in Figure 4. Each site (node) has 14 incident bonds (edges): eight normal to hexagonal 
faces (nearer neighbours) denoted by Bi with length Li = = '/Qa and six normal 

to square faces (further neighbours) denoted by B 2 with length L 2 = S = y/Sa. We keep 
the expressions with both S and a, because S is more convenient when the 3-complex is 
sized using an average grain volume, S = (2Vgr)'^ ■, while a is more convenient when it 
is sized using an average triple line length, a = Ltl, in the real microstructure. Note, 
that the construction of the graph follows the construction of a 3-complex dual to the one 
representing the microstructure, but it terminates before introducing 2-cells dual to the 
1-cells (edges) and 3-cells dual to the 0-cells (vertices) of the original 3-complex. 


2 The model with calculus on discrete mani¬ 
folds 

In the discrete calculus theory the boundary of a p-complex is represented by the so-called 
incidence matrix, see [5]. Specifically for 1-complexes (graphs), the incidence matrix, takes 
the simplest form denoted hy A = [aij]iZi’ 2 ’.'.'.’m- Where 

{ 0, if node j is not a node of edge i I 
1, if node j is the first node of edge i > . 

— 1, if node j is the second node of edge i ) 

For the the node-edge incidence matrix A G where m is the number of edges and 

n is the number of nodes of the graph. Row i, corresponding to edge i, contains exactly 
one 1 in the column corresponding to the first node of the edge and one -1 in the column 
corresponding to the second node of the edge. All remaining coefficients of the matrix 
are zeros. The incidence matrix contains not only the topology of the graph, but also 
represent the boundary operator. Thus application of A € to a nodal function is an 
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(Edge, 



B,(4) 


Figure 4: Illustration of site-bond model construction as a reduction of 3- 
complex to 1- complex (graph). Left figure shows a 3-cell with two neighbor¬ 
ing 3-cells which goes into 3 sites (0-cells, graph nodes) connected by bonds 
(1-cells, graph edges). Right hgure shows a single 3-cell with the 14 bonds 
(1-cells, graph edges) in the substituting site-bond topology. The two types 
of bonds, denoted by Bi and i? 2 , have distinct lengths Li and L 2 . 


edge function - the gradient of the nodal. Further, application of the transposed incidence 
matrix to an edge function is a nodal function - the divergence of the edge. Since we will 
know the initial positions of the nodes and the edges from a partiocular cellular structure, 
the incidence matrix will be known. 

Assume a lattice of n nodes connected through m edges. Let Di S i = l,2,...,n 
be a discrete (vector-valued) function over nodes and bi S i = l,2,...,m a discrete 
function over edges. The relation 

AD = b. 

states that b is the gradient of D, where 
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where Xi S z = 1, 2, n are the nodal coordinate and Ui e z = 1, 2, m the edge 
coordinates, the modulus of which provide edge lengths. Then 

= y. (1) 

In Q, the components of y are unknown, but some components of X are known from 
the essential boundary conditions required to complete the boundary value problem. The 
loading on our system will be applied only via essential boundary conditions, i.e. displace¬ 
ment control, where the coordinates of some nodes change from initial to prescribed final 
values. 

Solid materials accommodate strain from external loading by reversible (elastic) rear¬ 
rangement, giving rise to internal stresses, and by dissipating energy via slip (plasticity) 
or separation (surface generation). In discrete settings it has been shown that the stress 
between two regions is a traction vector always normal to their interface m- In the graph 
framework, this means that the stress is a vector (force, not a tensor as in continuum 
mechanics) acting along the current (deformed) orientation of an edge. Let rii = |^, 
z = 1,2, ...,77z, be the unit vectors along edges, where \yi\, i = 1,2, are the edge 

lengths after some deformation. The edge forces Fi, z = 1,2, ...,m, are therefore given by 

Fi = \Fi\ni=^-^yi, Vz = 1, 2,...,m. 

m 

This can be summarised by 

F = K{y)y, (2) 

where 

■ l^’il 
l^’2| 


l^ml 

Let \bi\, i = l,2,...,m, be the initial lengths of bonds. Then |Fi|, i = l,2,...,m, 
are related to the edge elongations, |6i| — \yi\, via potentially a non-smooth function as 
illustrated in Figure 5. 

Bond elongation, positive x-axis in Figure 5, results in positive or tensile force (or zero 
after x^i). Bond contraction, results in negative or compressive force. Bond elongation 
is related to the established in mechanics notions of stretch and strain. Stretch is the 
ratio between deformed and initial lengths, i.e. and strain is defined as the natural 
logarithm of stretch in the general case. For stretches around 1, the strain is approximated 
by the engineering strain . Mathematically the stretch can be anything from zero to 

one under contraction and from one to infinity under elongation. In practice, the possible 
stretches depend on the material in question. Rubber for example allows for stretches 
more than 2 (more than 100 percent elongation) without any plasticity or failure. Steels 
become plastic at strains between 0.2 percent and 0.5 percent. Hence the difference 
is really small at the onset of plasticity, represented by the line with zero slope in Figure 
5. Surface separation, represented by the line with negative slope (softening branch), is 


K{y) = diag \ —-,....-— 

|yi| ml Wm 


(3) 
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Figure 5: Conceptual edge forces-displacement relation. 


initiated at several percent strain and ultimate failure occurs typically between 40 percent 
and 60 percent. 

Let Bi, i = 1,2, ...,n be external forces at nodes, given as natural boundary condi¬ 
tions. These can all be zero as in the case of displacement-controlled loading, but can 
in general represent force-controlled loading. The balance of linear momentum at nodes 
(force equilibrium) requires that 

A^F = B. (4) 

where S j^nxm jg transpose of the incidence matrix A and 


B = 


Bi 

B2 

Br, 


Dnx3 


In summary, at each node we have a position vector (3 components in three coordinate 
directions) and a force vector (3 components in the same directions). The discrete position 
field on nodes is X, containing 3n numbers, where n is the number of nodes. The discrete 
force field on nodes B, contains 3n numbers. Applying AX gives us the discrete position 
field on edges y. This is a field containing 3 to numbers, three coordinate differences per 
edge. The position field on edges gives rise to a force field in the edges, which resists the 









deformation. Thus is a discrete force field on the edges with 3m components. Hence we 
have F = f{y) = K{y)y, which is a material-dependent function defining the resistance to 
deformation. Applying gives us a discrete force held on the nodes. This must equal 
B to balance the linear momentum in the system, i.e. ensure equilibrium. 

We note the exclusion principle for nodal helds in the system dehnition: an element of 
B is known (the force in a given direction at a given node) if and only if the corresponding 
element of X is unknown (the displacement in the same direction at the same node). This 
means that when a nodal force is prescribed (a natural boundary condition), the system 
solution provides its conjugate displacement. Inversely, when a nodal displacement is 
prescribed (essential boundary condition), the system solution provides its conjugate force. 

Substitution of ([^ into Q gives 

A^K{y)y = B, 

whereby using 0 into the above expression provides the formal description of the me¬ 
chanical deformation of the graph with 

A^K{AX)AX = B. (5) 


3 Solving the non-linear system 

System ([^ is a non-linear system. This is because of the diagonal matrix K defined in 
0. We provide the following Theorem 


Theorem 3.1. Consider the system ([^. Then 
(a) An effective linearization of system 0 is 

Ax = B. 


Where A = A^KA and for 0 < e << 1 


K = 




'-2xoi 


\b.\ 

+Xoi ' 

b.\ 


fii(x3i-\bi\) _ ^ 

^ tjl<2<m? 


diag[ 


2{x3i + \bi\){x2i + \bi\) 


Xoi < \yi\ - \bi\ < 0 
0 < |y*| - \bi\ < Xu 
Xu < \yi\ - \b,\ < X2t 
X2i < \yi\ - \h\ < X3i 


( 6 ) 


(b) Let 

^ _ ^11 Ai2 

A21 A22 

Then system 0 can be divided into the following subsystems 

AiiXp = Bp — Ai2Xq (7) 

and 

Bq = Al^Xp + A22Xq. ( 8 ) 

From the above systems only Q has to be solved. Then Xp can be replaced in (|^ 
and Bq is easily computed. 
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(c) For the solution of 0 

(i) If Ai 2 is full rank, then 

Xp = AY,^(Bp-Ai2Xg). (9) 

(ii) If All is rank deficient, then an optimal solution Xp for Q is given by 

Xp = (ifiiii + E^E)-^A^i{Bp - Ai2Xg). (10) 

Where E is a matrix such that A^iAn + is invertible and ||i ?||2 = 

0 < 0 << 1. Where 11-112 is the Euclidean norm and with we denote the 
transpose tensor. 

Proof. For the proof of (a) we consider the system (1^. Vi = 1,2, we will seek 
I F I 

bounds for We have the following cases. 


For xoi < \yi\ - |6*| < 0 we have |F,| = - Ihl), 

l-Fil _ _ \^i \ \ 

\yi\ xoi \yz\ 

Since ccoi < \yi\-\bi\ < 0, we have a;oi + |fei| < \yi\ < 

or, equivalently, 1 < |^ < or, equivalently, 


see Figure 5. Or, equivalently, 


( 11 ) 

or, equivalently, ^ 2 |±M < ||i| < 1, 


\h\ 

xoi + \bi\ 


h\ 

\yi\ 


< 0 . 


Hence, by using the above inequality and 0 we get 

|Fi| ^ foi _ \bi\ ^ ^ ^ 

\yi\ \bi\+xiii 


For 0 < \yi\ - \bi\ < Xn we have \Fi\ = ^^{\yi\ - \bi\), see Figure 5, or, equivalently. 


|Fi| _ /h/i _ N 
\y^\~Xl^ \y,\^- 

Since 0 < \yi\ — \bi\ < xu, we have 

0<l-f^<l- 

\yi\ xii + \bi\ 


Hence, by combining the above inequality with (12) we get 


\yi\ 


fli 

2xii 


(1 


\h\ 

\bi\ + Xu 


)±e. 


( 12 ) 
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For Xu < \yi\ — |6i| < X 2 i we have \Fi\ = fu, see Figure 1, or equivalently, 


\yi\ \yi\' 


( 13 ) 


Since xu < |j/i| - \bi\ < X21, we have xu + |&i| < \yi\ < X2i + |&i|, or equivalently 

1 ^ 1 ^ 1 

X21 + 1^*1 “ \yi\ ~ Xu + l^il 

or, equivalently, by using the above inequality and we get 

/h < ^ < /h 

X2i + 1^*1 “ \yi\ ~ Xu + \hi\' 


Hence 


For X 2 i <\yi\ 
equivalently. 


\Fi\ ^ fuix2i - Xu) ^ 
ly^l 2 {\bi\ + xu){\h\ + X2i) 

|5,| < X 3 , we have |F,| = - |6*|) + see Figure 5, or 

(i-N)+ A- ( 14 ) 

\yi\ X2i-X3i \yi\ X3i-X2i\yi\ 


Since X21 < \yi\ - \hi\ < xai, we have X2i + \bi\ < \yi\ < xst + \bi\, or equivalently, 
^ jfll ^ 01' equivalently 


h\ < N < 


or, equivalently. 

Similarly, 


X3i + \bi\ \yi\ X2i + \bi\' 


1 _ < 1 _ M < 1 _ JL. 

X2i + \bi\ ~ \yi\ ~ X3^ + \bi \■ 


1 1 1 
< ^-r < 


X3i + \bi\ \yi\ X 2 i + \bi\' 

Hence, by using the above two inequalities combined with (14) we have 

, \b^\ 




fli 


\bi\ ^ ^ X 3 ifu 


\yt\ 2{x2t - X3i) \bi\ + X 2 i \bi\+X3i 2{x3i - X 2 i) \bi\ + X 2 i \h\+X3i 


±e, 


or, equivalently, 

\Fi\ ^ hi \bi\ {X 3 i - 1 ) ^ ^ 

\yt\ 2{\bi\ + X2i){\b^\ + X2i) 

For the proof of (b) system (|^ can be written as 


All 

^12 


■ Xp ■ 


Bp 

A21 

A22 


. . 


. B, . 
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or, equivalently, 


+ Al2^q — Bp, 

-AlsXp + ^22^5 = Bq. 

From the above expressions we get the subsystems 0, (§. 

For the proof of (c), for (i), since An is full rank, we arrive easy at 0. For (ii), since 
All is rank deficient, if [Bp — Ai 2 Xq] ^ colspanAn system Q has no solutions and if 
[Bp- Ai 2 Xq] £ colspanAn system Q has infinite solutions. Let 


such that the linear system 
or, equivalently the system 


L{Xp) =L + EXp, 

AnXp = L(Xp), 
(An - E)Xp = L 


has a unique solution. Where E is a, matrix such that A^iAiiEE'^E is invertible, ||if ||2 = 
9, 0 < 0 « 1 and EXp is orthogonal to L — AnXp. Hence we want to solve the following 
optimization problem 


or, equivalently, 

or, equivalently. 


mm 


(Bp - Ai2Xq) - L 
s.t. (An — E)Xp = L, 


(Bp - Ai 2 Xq) - (in - E)Xp 


mm 


(Bp — Ai2Xq) — AnXp 


EXr 


The reason for using the matrix E is because the matrix An is rank deficient and hence 
the matrix iniii is singular and not invertible. To sum up, we seek a solution Xp 
minimizing the functional 


Hi(Xp) = 


Bp — Ai2Xq — AnXp 


+ 


EXr 


Expanding Hi(Xp) gives 

Hi(Xp) = (Bp - Ai2Xq - AnXpf(Bp - i^X, - AnXp) + (EXp)^EXp, 
or, equivalently, 

Hi(Xp) = (Bp-Ai2Xq)^(Bp-Ai2Xq)-2(Bp-Ai2Xq)'^AnXp+(Xp)'^AliAnXp+(Xp)'^E'^EXp 
because (Bp — Ai 2 Xq)'^AnXp = (Xp)'^An(Bp — A^Xq). Furthermore 

\T 


d 


Hi(Xp) = -2Aii(Bp - Ai2Xq) + 2AiiAnXp + 2E^ EXp. 
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Setting the derivative to zero, -^^Hi{Xp) = 0, we get 

(^n^ii + E)Xp = — Ai2Xq). 

The solution is then given by 

Ip = (ifiiii + E^E)-^Al^{Bp - ii2X,). 

Hence an optimal solution of 0 is given by ([To|). The proof is completed. 


4 Numerical example 



Figure 6: A unit cell 


We assume a unit cell of 15 nodes and 14 edges. Let Di, i = 1,2,..., 15, bi and \bi\, 
xo^, Xu, X2i, X3i, foi, fii, i = 1,2,..., 14, for * = 1,2,..., 14, be given by Table 1. Then 
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In addition, if X = 
be unknown and Xq, 


x 3 


B = 


Xp 
Xq 

Bp given as follows: 


Bp 

B, 


G with p 


7, q = 8, let Xp, 
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^ D'l bi \bi\ Xqi X\i X2i X^i foi fli 


1 

( 0 , 0 , 0 ) 

( o ,- i , i ) 


- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

2 

( 0 , 1 , 0 ) 

( 0 , 1 , 0 ) 

1 

- 0.2 

0.2 

0.3 

0.5 

- 0.1 

0.1 

3 

( 0 , 0 , 1 ) 

( 1 ,- 1 ,-!) 


- 0.1 

0.1 

0.2 

0.3 

- 0.2 

0.2 

4 

( 0 , 1 , 1 ) 

( 0 , 1 , 0 ) 

1 

- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

5 

( 1 , 0 , 0 ) 

( o ,- i , i ) 


- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

6 

( 1 , 1 , 0 ) 

( 0 , 1 , 0 ) 

1 

- 0.2 

0.2 

0.3 

0.5 

- 0.1 

0.1 

7 

( 1 , 0 , 1 ) 

( 1 ,- 1 ,- 1 ) 

73 

- 0.1 

0.1 

0.2 

0.3 

- 0.2 

0.2 

8 

( 1 , 1 , 1 ) 

( 0 , 1 , 0 ) 

1 

- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

9 

( 2 , 0 , 0 ) 

( 0 ,- l , l ) 

72 

- 0.2 

0.2 

0.3 

0.5 

- 0.1 

0.1 

10 

( 2 , 1 , 0 ) 

( 0 , 1,0 

1 

- 0.1 

0.1 

0.2 

0.3 

- 0.2 

0.2 

11 

( 2 , 0 , 1 ) 

( 0 , 1 , 0 ) 

1 

- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

12 

( 2 , 1 , 1 ) 

( 0 ,- l , l ) 

72 

- 0.2 

0.2 

0.3 

0.4 

- 0.2 

0.2 
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( 2 , 2 , 1 ) 

( 0 , 1 , 0 ) 

1 

- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

14 

( 2 , 1 , 2 ) 

(- 2 ,- l ,- 2 ) 

3 

- 0.1 

0.1 

0.3 

0.4 

- 0.1 

0.1 

15 

( 2 , 2 , 2 ) 










Table 1: Nodes, edges and other given values 


(a) For xq, < \y,\ - \bi\ < 0 


0.1 

1.1 

0.1 

0.1 

0.1 

1.1 

0.1 

1.1 

1.1 

1.1 

0.1 

0.1 

1.1 

1.1 

0.1 

1.1 

0.1 

1.1 

1.1 

1.1 

1.1 

1.1 

0.1 

0.1 


(b) For 0 < \yi\ - \h\ < xu 


0.2 

1.2 

0.2 

0.2 

0.2 

1.2 

0.2 

1.2 

1.2 

1.2 

0.2 

0.2 

1.2 

1.2 

0.2 

1.2 

0.2 

1.2 

1.2 

1.2 

1.2 

1.2 

0.2 

0.2 


0.1 

1.1 

0.1 

0.1 

0.1 

1.1 

0.1 

1.1 

1.1 

1.1 

0.1 

0.1 

1.1 

1.1 

0.1 

1.1 

0.1 

1.1 

1.1 

1.1 

1.1 


0.2 

1.2 

0.2 

0.2 

0.2 

1.2 

0.2 

1.2 

1.2 

1.2 

0.2 

0.2 

1.2 

1.2 

0.2 

1.2 

0.2 

1.2 

1.2 

1.2 

1.2 
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(c) For Xu < \yt\ - \hi\ < X 2 i 





0.3 

1.3 

0.3 

0.3 

0.3 

1.3 

0.3 

1.3 

1.3 

1.3 

0.3 

0.3 

1.3 

1.3 

0.3 

1.3 

0.3 

1.3 

1.3 

1.3 

1.3 

1.3 

0.3 

0.3 


0.3 

1.3 

0.3 

0.3 

0.3 

1.3 

0.3 

1.3 

1.3 

1.3 

0.3 

0.3 

1.3 

1.3 

0.3 

1.3 

0.3 

1.3 

1.3 

1.3 

1.3 


(d) For X 2 i < \yi\- \bi\ < x^i 




0.4 

1.4 

0.4 

0.4 

0.4 

1.4 

0.4 

1.4 

1.4 

1.4 

0.4 

0.4 

1.4 

1.4 

0.4 

1.4 

0.4 

1.4 

1.4 

1.4 

1.4 

1.4 

0.4 

0.4 


0.4 

1.4 

0.4 

0.4 

0.4 

1.4 

0.4 

1.4 

1.4 

1.4 

0.4 

0.4 

1.4 

1.4 

0.4 

1.4 

0.4 

1.4 

1.4 

1.4 

1.4 


From Theorem 3.1 the solution of systems Q, (|^, will be 
(a) For xqi < \y^\ - \hi\ < 0 




0 

111.2090 

111.0840 

107.6199 

101.1199 

84.1493 

79.7743 


0 

110.9590 

110.8340 

90.0494 

83.5494 

66.5788 

63.4538 


0 

152.1000 

150.7250 

112.6199 

101.1199 

84.1493 

79.7743 


,Bq 


-4.6000 

0.0000 

- 0.1000 

0.1000 

0 

0 

0 

0 


-3.4000 

-0.7657 

0.6657 

-0.3000 

0.4828 

-0.4828 

0.2667 

-0.0667 


-4.8000 

0.2000 

0.1000 

- 0.1000 

-0.2828 

0.2828 

0.0667 

-0.0667 


(b) For 0 < \yi\ - |6j| < xu 


Xp = 


0 

-129.8438 

-129.5938 

-122.6656 

-114.6656 

-94.8666 

-89.8666 


0 

-127.5938 

-127.3438 

-103.0951 

-95.0951 

-75.2961 

-71.5461 


0 

-170.7349 

-169.2349 

-127.6656 

-114.6656 

-94.8666 

-89.8666 


, Bn 


-5.2000 

0 

0.1000 

- 0.1000 

0.0000 

0 

0 

0 


-4.4000 

0.7657 

-0.6657 

0.3000 

-0.4828 

0.4828 

-0.2667 

0.0667 


-5.0000 

- 0.2000 

- 0.1000 

0.1000 

0.2828 

-0.2828 

-0.0667 

0.0667 
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(c) For Xu < \Vt\ - \bi\ < X 2 i 


0 

0 

0 

0.3591 

1.3559 

0.3749 

0.3581 

1.3549 

0.3708 

0.3564 

1.3504 

0.3634 

0.3431 

1.3371 

0.3431 

0.3308 

1.3248 

0.3308 

0.3164 

1.3136 

0.3164 


(d) For X 2 i < \yz\ - \bi\ < x^i 


-5.8 

138.2 

-148.8 

0 

-696.42 

143 

-132 

685.42 

132 

132 

-275 

-132 

0 

419.71 

-276.71 

0 

-419.71 

276.71 

0 

1166 

1023 

0 

-1023 

-1023 


0 

0 

0 

-0.026822 

0.986 

-0.12908 

-0.021694 

0.99113 

-0.11113 

0.0074869 

1.0568 

-0.0089966 

0.043751 

1.093 

0.043751 

0.10245 

1.1517 

0.10245 

0.16655 

1.203 

0.16655 


-6.4 

-66.067 

54.267 

0 

132.45 

-60.667 

22.286 

-94.071 

-22.286 

-22.286 

82.952 

22.286 

0 

-91.33 

30.664 

0 

91.33 

-30.664 

0 

-146.97 

-86.308 

0 

86.308 

86.308 


In addition, the edge elongations are given by Tables 2, 3. 


Vi 

xoi < \ yi \ - \ bi \ < 0 

0 < \ yi \ - \ bi \ < Xu 

Vi 

(0, 0, 0) 

(0, 0, 0) 

y2 

(-0.1250, -0.1250, -1.3750) 

(0.25, 0.25, 1.5) 

ys 

(-3.4641, -20.7846, -38.1051) 

(6.9282, 24.2487, 41.5692) 

yi 

(-6.5, -6.5, -11.5) 

(8 ,8, 13) 

y5 

(-16.9706, -16.9706, -16.9706) 

(19.799, 19.799, 19.799) 

y& 

(-4.375, -3.125, -4.375) 

(5, 3.75, 5) 

y7 

(-79.6743, -62.3538, -79.6743) 

(90.0666, 72.7461, 90.0666) 

ys 

1—1 

1—1 

o ' 

(0, -1, 1) 

1/9 

(0, 1, 0) 

(0, 1, 0) 

yio 

(1,-1, -1) 

(1, -1, -1) 

yii 

(0, 1, 0) 

(0, 1, 0) 

yi2 

(0, -1, 1) 

(0, -1, 1) 

yi3 

(0, 1, 0) 

(0, 1, 0) 

1/14 

(0, -1, -1) 

(0, -1, -1) 


Table 2: Edge elongations 
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Vi 

Xu < \yi\ - \hi\ < X 2 i 

X 2 i <\yi\- \bi\ < X 3 i 

yi 

(0, 0, 0) 

(0, 0, 0) 

y2 

(-0.00096154, -0.00096154, -0.0041667) 

(0.0051282, 0.0051282, 0.017949) 

ya 

(-0.0016951, -0.0045203, -0.0073454) 

(0.029181, 0.065657, 0.10213) 

yA 

(-0.013287, -0.013287, -0.02028) 

(0.036264, 0.036264, 0.052747) 

yb 

(-0.012328, -0.012328, -0.012328) 

(0.058701, 0.058701, 0.058701) 

y& 

(-0.014423, -0.011218, -0.014423) 

(0.064103, 0.051282, 0.064103) 

yr 

(-0.016386, -0.013561, -0.016386) 

(0.23345, 0.19697, 0.23345) 

ys 

(0,-1, 1) 

(0,-1, 1) 

y9 

(0, 1, 0) 

(0, 1, 0) 

yio 

(1,-1, -1) 

(1, -1, -1) 

yn 

(0, 1, 0) 

(0, 1, 0) 

yi2 

(0, -1, 1) 

(0, -1, 1) 

1/13 

(0, 1, 0) 

(0, 1, 0) 

yii 

1—1 

1—1 

o' 

(o,-i,-i) 


Table 3: Edge elongations 


Conclusions 

In this article we focused on the mathematical derivation of the behavior of lattice el¬ 
ements, nodes and edges, in the presence of non-linear deformations. We presented an 
efficient integrated lattice formulation and demonstrated the methodology with a simple 
example. An extension of this work is planned to test and verify the method, by com¬ 
paring simulations of lattice elastic-plastic-damage behavior to experimentally measured 
material responses. 
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